Quantum Monte Carlo for minimum energy structures 
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We present an efficient method to find minimum energy structures using energy estimates from 
accurate quantum Monte Carlo calculations. This method involves a stochastic process formed from 
the stochastic energy estimates from Monte Carlo that can be averaged to find precise structural 
minima while using inexpensive calculations with moderate statistical uncertainty. We demonstrate 
the applicability of the algorithm by minimizing the energy of the H2O-OH - complex and showing 
that the structural minima from quantum Monte Carlo calculations affect the qualitative behavior 
of the potential energy surface substantially. 



First principles electronic structure methods have been 
used to describe and explain a wide range of properties 
for different condensed matter systems. A critical step 
is the accurate determination of the ground state atomic 
structure, since many important properties of a mate- 
rial can change dramatically depending on the structure. 
Because of the balance between accuracy and computa- 
tional cost, density functional theory (DFT) has become 
a commonly used method to find equilibrium geometries 
of both molecules and extended systems. The primary 
reason for this is the availablility of forces with little ex- 
tra computational cost over the energy calculation. Using 
a typical quasi-Newton minimization algorithm, the local 
minimum of a potential energy surface can be found in 
O(Ndof), where Ndof is the number of degrees of free- 
dom to be optimized. This favorable scaling has made 
it possible to find minima for many systems of inter- 
est. However, in many situations, including transition 
metals, excited states, and weak binding, current density 
functional theories may not be accurate enough even for 
structures, and more accurate post-Hartree-Fock meth- 
ods that scale as 0(N^~ 7 ), where N e is the number of 
electrons, can often be too computationally expensive. 

Quantum Monte Carlo (QMC), a stochastic approach 
to solving the many-body Schrodinger equation, offers 
favorable scaling to obtain the total energy, 0(N^~ 3 ), 
and has been shown to provide near chemical accuracy 
in many different systems ^H]- However, there are two 
major challenges in using QMC methods to obtain high 
precision minimum energy structures. The first is that 
the techniques so far proposed to calculate forces in diffu- 
sion Monte Carlo all have large variances and error terms 
that depend on the quality of the trial wave function, 
which is often poor in systems where DFT fails and one 
would like to apply QMC methods. In fact, despite much 
work in recent years [SHE], QMC forces using the highly 
accurate diffusion Monte Carlo method have not been 
applied to more than a few degrees of freedom, although 
the simpler variational Monte Carlo technique has been 
applied to more [9j . The second challenge is the stochastic 
nature of Monte Carlo algorithms, which provides uncer- 
tainty in any estimator that only decreases as the square 



root of the computer time. Reducing the uncertainty 
enough to resolve the minimum structure accurately us- 
ing forces or total energies is often prohibitively expensive 
computationally. Methods such as the stochastic gradi- 
ent approximation [TUl [TT] that are able to operate in the 
presence of noise suffer from this large uncertainty in the 
forces. As a result, there are no geometry optimizations 
of more than three [12] degrees of freedom, to our knowl- 
edge. 

In this article, we describe an algorithm that uses the 
already-accurate total energies from QMC to obtain min- 
imum energy geometries with well-defined stochastic un- 
certainties with multiple degrees of freedom. The algo- 
rithm consists of two major parts. One is a sequence of 
minimizations along one dimensional lines. The use of ID 
minimizations allows us to use efficient fits to determine 
the minimum precisely. The second part is a quadratic 
fit of the many-dimensional energy surface to determine 
the new search directions. Both of these parts are com- 
pletely aware of the stochastic uncertainty present in the 
evaluations of the energy, an important feature obtained 
by the use of Bayesian inference. We apply this approach 
to the hydrogen-transfer model of H 2 0-OH~ and show 
that our method can help clarify challenging problems 
that require accurate calculation of the electronic ground 
state. 

DFT and Hartree-Fock calculations were performed 
using the GAMESS[I3] package. We used soft 
pseudopotentials[14] to remove the core electrons and a 
TZP-quality contracted gaussian basis to represent the 
one-particle orbitals. All-electron calculations were per- 
formed using the aug-cc-pVQZjTS] basis to check the ba- 
sis set and pseudopotential errors. All QMC calculations 
were performed using the QWalk[16 program. For en- 
ergies, we used fixed node diffusion Monte Carlo (DMC) 
with a time step of 0.02-0.05 Hartrees -1 , converged for 
the properties of interest. The trial function was a Slater- 
Jastrow function with hybrid PBE0[T7] orbitals, one and 
two body terms in the Jastrow, and further checks of 
the localization error with a three-body Jastrow factor. 
Further details can be found in e.g. Refs [IB] US]. 

Our minimization algorithm is similar in spirit to many 
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FIG. 1: (color online) An outline of the algorithm. See text 
for details. 



successful minimization algorithms in that it is based on 
minimizing along directions and updating an estimate of 
the Hessian matrix to find the diagonal directions. How- 
ever, in our approach the Hessian matrix is inferred using 
the Bayesian interpretation of probability, which has the 
effect of making the algorithm very robust to stochas- 
tic uncertainty. Here we present the case when only the 
value can be calculated, but gradients of the objective 
function can be included easily if they are available, in- 
creasing the efficiency. By using inference, we are able to 
make very efficient use of the available data to find the 
Hessian matrix without having to reduce the stochastic 
uncertainties to small values that cost large amounts of 
computational time to achieve. 

The kernel of the algorithm is a sequence of line min- 
imizations. We show that a sequence of uncertain line 
minimizations will obtain the true minimum on aver- 
age. Such a sequence can then be viewed as a gen- 
erator of random variables whose expectation value is 
the true minimum. Suppose on the first step, we start 
with the approximate minimum x' ' . We then define 
(5x(°' = x^ ' — m, where m is the vector of the unknown 
true minima. We wish to design our process such that 



the expectation value of Sx^ equals zero as n — >• oo. 
Within the quadratic region around the minimum, the 
potential energy surface is given by 



E(x) = 5x T H5x + E a , 



(1) 



where H is the symmetric Hessian matrix and E is the 
minimum total energy. On minimizing along each direc- 
tion i, there are two components of the distance from 
the true minimum. The first is deterministic and comes 
about from non-zero Sxj for j ^ i. The second is the 
stochastic error from the uncertainty in the line mini- 
mization, which we can estimate using the Bayesian tech- 
niques above. We find that 
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where x\ 1S a random number. For the nth iteration, 
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The minimum along each line is found using line fitting, 
as described in the EPAPS document. This allows for a 
very efficient determination of the minimum. 

One can see from Eqn [3] that the smaller the off- 
diagonal matrix elements of the Hessian are, the less in- 
terference directions have on each other-for a diagonal 
Hessian, only one minimization in each direction is neces- 
sary once we are in the quadratic regime. We can use the 
information from the line minimizations to estimate the 
Hessian as the minimization proceeds. We parameterize 
the quadratic region with a set of parameters cq , includ- 
ing the elements of the Hessian matrix, the minima, and 
the minimum energy. After performing the line fits, we 
have distributions of line fit parameters, each given by C£ 
for line t. Using Bayes' theorem, the likelihood function 
of the quadratic parameters L is given by 

L[cq\D)<xp{D\cq) = H [ p(Dt\L t )p(Li\cQ)dci, (4) 

where D is the set of function evaluations and stochas- 
tic uncertainties given by e.g. QMC and each Di is the 
subset of function evaluations used to minimize along a 
given line I. Since p(Li\cq) is not based on stochastic 
data, it is a delta function that forces eg to be consis- 
tent with cq as follows. Since in the quadratic region, 
E(x) = Sx. T HSx + Eq, minimizing along a direction v 
from a starting position x gives the one-dimensional 
function of t, the position along the line: 

E(t) = (x +tv- m) T H(x + tv - m) + E . (5) 

This gives the following constraints on each set of line 
parameters for the minimum, curvature, and minimum 
function value cf. 
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The net result of this transformation is a set of parame- 
ters that includes the Hessian matrix, the location of the 
minimum, the objective function at the minimum, and 
the parameters used for the line fits above the three con- 
straints implied by the cq's. We examine the properties 
of this probability distribution function in the EPAPS 
material. It turns out that the maximum likelihood esti- 
mator is typically accurate enough to determine the Hes- 
sian, which is then diagonalized to obtain new search 
directions. 

We outline a single iteration of the algorithm in Fig [T] 
This step operates in two principle regimes. The first is 
far away from the minimum. In this regime, the Hessian 
inference method behaves similarly to a deterministic di- 
rection set method, with the directions being determined 
by the Hessian inference. The second regime is when the 
calculation has mostly converged and the deterministic 
error in Eqn[2]is small compared to the stochastic error. 
In this regime, deterministic direction choices are partic- 
ularly useless, but the inferred Hessian method is able 
to automatically account for the stochastic error. Once 
in this "stochastic regime," we can use Eqn [3] to justify 
averaging the x> n ' 's to obtain a more precise estimate of 
the minimum. For stochastic functions, the performance 
of the Hessian inference is much higher than traditional 
methods such as Powell's method[TS], which we compare 
in Fig [2] Since Powell's method uses the concept of con- 
jugate directions to find the search directions, the un- 
certain minimizations cause the algorithm to fail quickly. 
Even when the directions are reset every sweep, they can 
become linearly dependent within a single sweep, causing 
a high failure rate for even moderate dimensionality. 

To show the value of optimizing geometries within 
QMC, we apply the method to the H20-OH~ complex 
(Figp3J), which is present in liquid water and important in 
many systems in condensed matter, biology, and chem- 
istry. The shape of this potential energy surface is cru- 
cial to understanding hydrogen transfer in water. In this 
case, we use our knowledge of the system to choose fixed 
search directions for efficiency reasons, omitting the Hes- 
sian inference method. For systems that are not as easily 
decomposed, the Hessian inference is invaluable. It has 
been noted [20 that current DFT functional disagree on 
the ground state structure of this complex. The potential 
minima are the non-centrosymmetric structure (A) and 
the centrosymmetric structure (B). 

Hartree-Fock and second order M0ller-Plesset pertur- 
bation theory (MP2) find that structure A is lower in 
energy, with a barrier to transfer the proton. This is the 
traditional picture of this structure. The local density 
approximation and generalized gradient approximation 
(PBE [21]) of DFT find that structure B is lower in en- 
ergy. Using our QMC line minimization method without 
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FIG. 2: (color online) Estimates of reliability for the Hessian 
inference method versus Powell line minimization. There are 
several runs for each number of degrees of freedom, each with 
a different potential energy surface and relative starting po- 
sition. The potential surface was a randomly generated pos- 
itive definite quadratic surface with a random minimum and 
a small cubic perturbation. The condition number of the ma- 
trices was typically 50 or less. We added a Gaussian random 
variable to the values of the potential energy surface to simu- 
late uncertain evaluation. Convergence was judged when the 
RMS distance to the true minimum converged to the stochas- 
tic floor of the minimization. 



TABLE I: H 2 0'-OH". The non-QMC methods are optimized 
using conjugate gradient routines in GAMESS. The asterisked 
DMC result corresponds to constraining the DMC geometry 
search to the symmetry of structure B. 

Method O-O' O'-H Structure type 



LDA 


2.448 


1.224 
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PBE 


2.470 


1.235 
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MP2 


2.469 


1.123 
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DMC 


2.491(2) 


1.111(3) 
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DMC* 


2.469(3) 


1.235(2) 


B 



constraints (seven degrees of freedom) we find structure 
A to be the minimum energy. The results are summa- 
rized in Table [I] DMC differs qualitatively from the DFT 
results in that structure A is the minimum, and quanti- 
tatively with MP2, since the oxygen-oxygen distance in 
MP2 is much smaller than in DMC for structure A. 

In Fig [3| we present the DMC energies of several min- 
ima. Without the geometry optimization algorithm, we 
would use minima from some other methods, for example, 
PBE and MP2 to obtain structures B and A respectively, 
then evaluate the total energy using DMC. In this case, 
since the MP2 approximation obtains a poor 0-0 dis- 
tance, DMC predicts a much higher energy for structure 
A and thus favors the centrosymmetric geometry. How- 
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FIG. 3: (color online) The relative DMC energies of geome- 
tries obtained by minimizing different potential energy sur- 
faces. The line is a guide to the eye. Stochastic error bars are 
the size of the symbols. 



ever, this is qualitatively incorrect, as we can see from 
the DMC-optimized structures, which predict structure 
B to be about 0.015 eV higher in energy than structure 
A. It is notable that 0.015 eV is a very small energy dif- 
ference on the scale of chemical bonding, which is why a 
small error in the geometry by using a lower level theory 
is enough to reverse the ordering. This energy difference 
is an upper bound to the barrier to transfer the hydro- 
gen, and thus at room temperature the hydrogen is free 
to transfer between the oxygen atoms. This may have 
important implications for the development of effective 
solvent models for water. 

In this work, we have viewed the DMC potential en- 
ergy surface as a given without attempting to improve 
the accuracy over a simple Slater- Jastrow form; how- 
ever, our method can also be used to optimize the nodal 



surface directly and thus further improve the accuracy 
of the DMC calculation. Similarly, any Monte Carlo 
method that depends on continuous parameters could 
put the stochastic line minimization algorithm to use. 
Stochastic line minimization may find application in ex- 
periments whose objective is to obtain a precise minimum 
in many-dimensional space, but the experiments are diffi- 
cult to perform precisely. If it is instead easier to perform 
many iterations of the experiment, then the minimization 
method that we have outlined may be applicable. The 
scaling for a quadratic potential is 0(N^ )0F ) if only val- 
ues are used; if even approximate gradients are available, 
the scaling can be reduced to approximately O(Ndof), 
as discussed in the EPAPS document. 

The algorithm is quite general, requiring only evalua- 
tions of the objective function with a known distribution, 
and can thus be applied to many minimization problems. 
At its core, it is a method for converting noisy value cal- 
culations into precise minima with rigorous error bounds, 
so it can be applied in any of the many problems where 
that is necessary. In electronic structure, diffusion Monte 
Carlo is accurate and applicable to many systems beyond 
the ground state of molecules, including solid structures 
and excited states [35]. The scaling for DMC geometry 
optimization is then 0{N c ^q F N^~ 3 ), depending on the 
size of the system and whether gradients are available. 
This favorable scaling combined with the already known 
high accuracy of quantum Monte Carlo could open up 
new levels of accuracy in structure determination in con- 
densed matter and chemical physics. 

This work was supported by the Department of En- 
ergy under grant No. de-sc0002623 and the NSF under 
grant No. 0425914. We wish to thank NERSC and Ter- 
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discussions. 
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Appendix A: Details of the minimization algorithm 

In this EPAPS document, we step through the mini- 
mization process in detail. Note that some of the details 
are not necessarily perfectly optimal, but they do work 
robustly and efficiently. We thus provide this description 
only as a reference implementation that could potentially 
be improved significantly. 



1. Line minimization 

Suppose we wish to minimize a function along a line 
using the values. For a deterministic function, a good 
method is the golden ratio method. However, for a non- 
deterministic function, locating the minimum by compar- 
ing the value is bound to be inefficient, since the objective 
function f(x) scales as f(x) ex x 2 , and the work to re- 
solve energy differences scales as the inverse square root 
of the difference. Therefore, near the minimum, the ef- 
fort to resolve position differences scales as the inverse 
quartic root of the difference. The computational cost 
using fitting (and gradients) on the other hand, scales as 
the inverse square root. Whether fitting using the value 
of the function or the gradient of the function is more 
accurate will depend on their relative variances and the 
specific problem. In this work, we use exclusively values 
since they are easily attainable in quantum Monte Carlo 
methods. 

We developed a reliable strategy to bracket the mini- 
mum. Given a starting position (to) and a trust radius r t , 



the algorithm works with three points, initialized as fol- 
lows: ti = t ~r t , t 2 = t , t 3 = t + r t . Let fa = f(U)+\ 
be the function evaluated at ti with a normal random er- 
ror x with standard deviation a , which can be controlled 
by performing longer calculations. Define S(fi,fj) as 
being true when fi is significantly different from fj, i.e., 
about four standard deviations away. A basic algorithm 
as follows: 

1. Evaluate /i, f 2 , and f 3 at uncertainty a 

2. IfS(/i,/ 2 )andSCf 2 ,/ 3 ) : 

• If f 2 < h and f 2 < fi stop. We have brack- 
eted the minimum. 

• If h < h and h > A, then set t t = t t - r t 
for i = 1,2, 3, return to 1. 

• K h > A and f 2 < /1, then set t i = U + n 
for i = 1,2, 3, return to 1. 

3. Else a = a/2 and return to 1. 

This algorithm occasionally gets stuck when two points 
are roughly the same distance from the minimum and 
thus have nearly the same value. This can be amelio- 
rated by searching between points when only two are 
insignificantly different. 

Once the minimum is bracketed to within the trust 
region, we can then fill in points along the line and fit. 
The distribution of points is not very critical, so long as 
there are a few points far away from the minimum to 
determine the curvature and higher order parameters ac- 
curately. We typically use a four parameter cubic form: 
f(x) — c + c 2 (x — C1) 2 + c 3 (x — C1) 3 , which gives a good 
balance between being accurate, easy to fit, and few pa- 
rameters. A quadratic form requires trust radii that are 
very small, reducing the efficiency of the algorithm. We 
perform the fit using Bayesian inference given the set of 
data D = {ti, fi, <7i}. The likelihood function is given by 



L(c\D)<xp(D\c), (Al) 



where 



p(D\c) ex []exp[-(/(i,,c) - fi) 2 /2a 2 ] (A2) 

i 

is the probability that we would have obtained the data 
given a set of parameters. The distribution of a given 
parameter a can be found by taking the marginal dis- 
tribution of L(c\D), or for more speed at the cost of not 
knowing the uncertainty well, simply using the maximum 
likelihood estimator. 



2. Hessian estimation 

One begins by noting that near the minimum, the en- 
ergy is described by a quadratic function: 



E(x) = Sx T H5x + Eo 



(A3) 



Far away from the minimum, this is an approximation 
to speed up the minimum search. Near the minimum, 
we used this property to prove the average properties of 
the stochastic sequence of line minimizations. Recall that 
only the minimum of the line minimization must be in the 
quadratic region. With no further assumptions, we can 
estimate the Hessian (and full quadratic surface) given a 
collection of line minimizations, even in the presence of 
stochastic noise. 

Let cq be the collection of all parameters for the 
quadratic region, eg the collection of line parameters for 
line £. We wish to find p(cq\D), where D is the set of 
function evaluations, arranged in lines. By Bayes' theo- 
rem, we know that 

p(c Q \D) ocp(£>|c Q ) = [] f p(Dg\L e )p(Lg\c Q )dcg (A4) 

Since p{Lp\cq) is not based on stochastic data, it is a 
delta function that forces eg to be consistent with Cq as 
follows. In the quadratic region, E(x) = 5x T HSx + Eq, 
minimizing along a direction v from a starting position 
xo gives the one-dimensional function 

E(t) = (x +tv- m) T H{x + tv - m) + E . (A5) 

This gives the following constraints on each set of line 
parameters eg: 

v T H (x - m) + (x - m) T Hv 

t m = TfTTZ (A6) 



2v T ffv 



d 2 E 
dt 2 



2v T iJv 



(A7) 



E(t m ) = (x + t m v - m) T H(x + t m v - m) (A8) 

We thus only need to integrate over the variables in each 
line that are not specified by the above equations. We can 
replace the eg by the reduced line parameters bg. Since 
there are three constraints per line, the dimensionality 
of eg is three less than the dimensionality of eg. We can 
recover the full set of line parameters as a function of the 
quadratic parameters and the reduced line parameters; 
i.e., Cg(c Ql c e ) 

We eventually want the likelihood function for the 
quadratic parameters and the reduced line parameters: 
L(cq, {cg}\D), where D is the data. This is given by 



L(c Q ,{c e }\{Dg}) ex Y[L{cg(c Q ,cg)\Dg). 



(A9) 



where L(cg\Dg) is given by Eqn Al 



The first estimate of the Hessian matrix is achieved 
by minimizing the violation of Eqn | A7| from the line fits. 
The minimum of the quadratic model is estimated as 
the last minimum found in a line minimization, and the 



minimum energy similarly. The parameters are then op- 
timized for maximum likelihood using a conjugate gra- 
dient optimization method. Finally, the parameters are 
randomized and reoptimized to get out of local maxima. 
We find that the maximum likelihood estimator is typi- 
cally quite sufficient for determining search directions; it 
is not necessary to sample to find distributions. 

Note that each line minimization places three con- 
straints on cq and cq has a dimensionality of 1+Noof + 
Ndof{Ndof + l)/2, so to make this a fully determined 
problem, 1/6 + N DO f/3 + N D qf{Ndof + 1)/6 line min- 
imizations are needed. This is always smaller than the 
Np OF line minimizations necessary in Powell's method, 
the standard in gradient-free minimization, with the ad- 
ditional benefit that the Hessian inference method is 
uncertainty-aware, which allows much higher efficiency 
when uncertainty can be traded for computational time. 

If gradient information is available, then they can 
be integrated directly into the likelihood function from 
Eqn |A9| as a simple product: 



L(c Q ,{cg}\{Dg},{gi}) oc 



Y[L(cg(c Q ,cg)\Dg) 



(A10) 



in* 



Q\&i 



where the likelihood function can be approximated as 
i(c Q |g,)cx (All) 

exph(V x £(x 4 , Cq) - g. i ) T a i (V XJ B(x 4 , Cq) - &)], 

and Ci is the covariance matrix between all the gradi- 
ent components. The Hessian inference then becomes 
essentially a quasi-Newton algorithm that is uncertainty 
aware. Approximate gradients are also useful, since min- 
ima can be determined by using only the value to perform 
the line fits. In some instances, one can imagine elabora- 
tions of this where DFT or other approximations to the 
gradient are used to accelerate the minimum search. 



3. A detailed look as the algorithm proceeds 

1. Start a line minimization sweep, for each direction 

• Bracket the minimum 

• Fill in data points within the trust region from 
the minimum 

• Fit the data points with a function 

• Update the minimum 

• Add fit to list of fits 

2. Prune list of fits to contain only the last Nhistory 
sweeps 

3. Fit N-dimensional quadratic model to lines, find 
maximum likelihood 

4. Diagonalize Hessian 




10 20 

Sweep number 



FIG. 4: Typical convergence of the eigenvalues and the 
estimate of the minimum for a nine-dimensional optimiza- 
tion problem. The potential energy surface was a random 
quadratic function with a cubic perturbation. The potential 
energy surface was calculated as E(x) + \i where x was a 
random Gaussian variable with mean zero. 



5. Update directions 

6. Return to 1. 

In Fig |4j we show the convergence of the eigenvalues 
from the Hessian estimation to the exact ones, and the 
minima from the line minimizations. 



Suppose that Ylj& \Hij/Hu\ < £,Vi for some < x < 1. 



Then|E fe ^^if|< 



H,, 



and so on. 



1. The average is equal to the true minimum 

Suppose that (x| j) ) = OVi, j. Then 

{5x[ n) )<x n 5xf. (B3) 

This goes to zero exponentially in n, so we can obtain an 
average arbitrarily close to zero by starting the average 
at high enough n. 

2. The average has finite variance 



Assume Var(vi) < u m Vi. Then 



Var(faf } ) < v m (l + x 2 + x 4 + ...) 



for large n. 



l-x 2 



(B4) 



3. The sequence has a finite correlation time 



Suppose that 



Appendix B: Proofs of the properties of the 
stochastic sequence 

We have the following equations: 

Jf-M _ v («) _ V^ H ij (n-l) 






37* 

HijHjk __ (n-2) 

HjiHi 



'Xk 



Hi 



Var(fo^) = V ar(xJ n) ) + E »l Var (*i ) 



&i 



-E v -(xr 2) ) 



(Bl) 



(xi^H^Wi,* 



(B5) 



for some variance a 2 . Then the correlation between ele- 
ments n and j in the sequence is 

\(6x\ ]) Sx[ n) )\ < a 2 E ^" _fc - tV"* E **' ( B6 ) 



fe=0 



fc=0 



where we reversed the sum indices. As n and j become 
large while keeping the difference finite, we can then find 
that 



Ea -<%)(!-%■) 



(B2)' 



\(Sx^Sx^)\<^-^ 



(B7) 



This goes to zero exponentially with n — j, and since the 
variance is finite, the correlation coefficient also goes to 
zero exponentially. 



